library(AUCell)
library(Seurat)
library(tidyseurat)
library(tidyverse)

dir.create("AUCell_tutorial")
setwd("AUCell_tutorial")

# read a sparse mtx
sobj_kera <-
  read_rds('../mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

exprMatrix <- sobj_kera |>
  GetAssayData()

exprMatrix |> glimpse()

# prepare gene sets -----------
# ?AUCell_calcAUC
## from GSEAbase --------
library(GSEABase)
genes <- c("gene1", "gene2", "gene3")
geneSets <- GeneSet(genes, setName="geneSet1")
geneSets

## character list ---------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

# score gene signatures -----------
# 21s for 12k cells
# 7s for 5 cores
# 4s for 78 cores
# use top 5% genes by default, increase maxrank if filtered data
system.time(
  cells_AUC <- AUCell_run(exprMatrix, kegg_mva, aucMaxRank = 2000,
                          BPPARAM=BiocParallel::MulticoreParam())
)

set.seed(333)
par(mfrow=c(1,1)) 
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, assign=TRUE) 

# thick vertical line: threshold for assignment
# selected from multiple thresholds
cells_assignment$geneSet$aucThr

selectedThresholds <- getThresholdSelected(cells_assignment)

# assigned "active" cells
mva.assigned <- cells_assignment$geneSet$assignment

mva.assigned |> glimpse()

sobj_kera |>
  mutate(mva.aucell = .cell %in% mva.assigned) |>
  ggplot(aes(trpv3_status, fill = mva.aucell)) +
  geom_bar(position = 'fill')

# plot result on tsne ----------
selectedThresholds <- getThresholdSelected(cells_assignment)

sobj_kera <- cells_AUC |>
  getAUC() |>
  as_tibble() |>
  pivot_longer(everything(), names_to = '.cell', values_to = 'aucell.mva') |>
  left_join(x = sobj_kera, y = _)

sobj_kera |>
  FeaturePlot('aucell.mva')

sobj_kera |>
  VlnPlot(c(kegg_mva, 'aucell.mva'),group.by = 'trpv3_status', pt.size = 0)

sobj_kera |>
  FindMarkers(group.by = 'trpv3_status', features = 'aucell.mva',
              ident.1 = 'Trpv3-high')

sobj_kera |>
  ggplot(aes(trpv3_status, aucell.mva)) +
  geom_boxplot() +
  ggpubr::stat_compare_means()

# compare with Seurat::AddModuleScore() ====
sobj_kera |>
  AddModuleScore(list(mva = kegg_mva), name = 'module.mva') |>
  ggplot(aes(trpv3_status, module.mva1)) +
  geom_violin(draw_quantiles = .5, scale = 'width') +
  ggpubr::stat_compare_means()

sobj_kera |>
  select(trpv3_status, aucell.mva) |>
  pivot_wider(names_from = trpv3_status, values_from = aucell.mva,
              values_fn = list) |>
  rowwise() |>
  summarise(pval = ks.test(x = `Trpv3-low`, y = `Trpv3-high`)$p.value)
